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Collective diffusion coefficient in a one dimensional lattice gas adsorbate is calculated using vari- 
ational approach. Particles interact via either a long-range, or a long range electron-gas-mediated 
(for a metallic substrate), or a 12 — 6 Lennard- Jones interaction. Diffusion coefficient as a func- 
tion of the adsorbate density strongly depends on the relationship between the substrate lattice 
constant and the characteristic length of the inter-particle interaction potential (which determines 
positions of the potential energy minima). The diffusion coefficient at fixed density as a function 
of the interaction characteristic length has an oscillating character due to the interplay between 
the inter-particle distances allowed by the substrate lattice structure and the average inter-particle 
distances which minimize the total interaction energy. 
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I. INTRODUCTION 

Manipulation of single atoms and self-assembly tech- 
niques are very useful methods in quantum engineer- 
ing. Self-assembly is determined by inter-atomic in- 
teractions, which at the crystal surfaces can be direct 
or indirect (i.e. induced by the substrate). Recent 
studies on adatom arrangement and their dynamics on 
metallic surfaces show that they experience an indirect 
electron-gas-mediated interactions. Electronic surface 
states are a source of a long range interactions, which 
decay with the inter-particle distance r as 1/r 2 and of- 
ten have an oscillatory character^^^i 6 -! 7 -!^. Such in- 
teractions between adatoms lead to their self-alignment 
in row a 8 ' 9 ' 10 ! 11 ; 12 ! 13 ) 14 or in hexagonal structures ^ 15 ! 16 . 
This ordering mechanism is a good candidate to be used 
for constructing and manipulating nanostructure sys- 
tems. Linear arrangements have unique magnetic and/or 
electronic properties. The ordering dynamics and a sta- 
bility of an ordered structure depend on diffusion of 
adatoms on the surface while diffusion in a system of 
adatoms is controlled by interactions between them. 

The collective or chemical diffusion coefficient of ad- 
sorbed species characterizes a relaxation of the local den- 
sity fluctuations in a many particle system. It involves 
jumps of individual atoms (adsorbed particles) from one 
binding site to another. Theoretical description of the 
collective diffusion process is a complicated many-body 
problem and various approaches have been applied to 
it, ranging from analytic ones based on master, Fokker- 
Planck, or Kramers equations to numerical Monte Carlo 
or molecular dynamics simulations. An important back- 
ground is provided in the works of Reed and Ehrlichia, an 
early summary by Comer—, and in reviews by Danani et 
al~2- and Ala-Nissila et ai 20 . The variational approach 
to the collective diffusion problem, used here, was de- 
veloped in a series of work a 21 ' 22 i 23 ' 24 i 25 ' 26 i 27 ' 28 and was 
shown to be a very efficient tool to analyze collective 
diffusion problems for various types of inter-particle in- 
teractions for homogeneous or inhomogeneous substrates 
either in one or two dimensions. 



We have shown in Ref. [28| that the long range re- 
pulsive interactions can be responsible for rapid macro- 
scopic rearrangements of adatoms upon minuscule adsor- 
bate density changes. We have shown that the adsorbate 
density dependent diffusion coefficient has peaks at den- 
sities corresponding to any ordered phase in a devil stair- 
case phase diagram. It is well know n 29 ' 30 that the devil 
staircase structure emerges when the inter-particle inter- 
action potential is repulsive and decays faster than 1/r. 
It has been shown, one the other hand, that adatoms on 
metallic surfaces often interact via forces with oscillat- 
ing in r potential energiesi^i 3 -^^^. They exhibit a 
1/r 2 decay modified by Friedel oscillations of the electron 
gas correlation function meaning that particles attract 
each other at some distances and repel at others. Such 
a distance dependence of the inter-particle interactions 
should, in principle, result in a diffusion character similar 
to that already investigated for pure 1 jr 2 repulsion. On 
the other hand, the oscillating interaction potential has 
local minima. Similarly, a minimum is present also for 
the interaction of the Lennard-Jones type. The question 
arises how the presence of such a minimum (or minima) 
affects the diffusion kinetics in a many particle system. 
In particular, it is an interesting question how does the 
oscillating character of the interactions decaying like 1/r 2 
modify diffusion in comparison with the interactions de- 
caying like 1 jr 2 monotonically. 

It is well known that the shape of interaction poten- 
tial, in particular, existence of its attractive parts is very 
important for the static behavior of the system. Sys- 
tems with attractive interactions form stable clusters of 
ordered phases, whereas ordering via repulsion is always 
global - it affects the entire system. Consequently, an 
ordered phase due to attractive interactions occurs for 
a wider range of adatom densities so it is easier to ob- 
serve experimentally 15 . Such difference in static proper- 
ties has to affect the dynamic behavior too, so it should 
affect also the diffusion process. For example, a fast col- 
lective diffusion in ordered structure a 27 ' 28 leads to a fast 
reorganization of the adsorbed layer. 

In this work we compare three types of interactions: 
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monotonically decaying like 1/r 2 , the oscillating ones de- 
caying like 1/r 2 with electron-gas-mediated oscillations, 
and the 12 — 6 Lennard- Jones interactions. We analyze 
the influence of the shape of the interaction potential on 
the diffusion coefficient at different adsorbate densities. 
Magnitude of the diffusion coefficient depends on sev- 
eral parameters, most of them related to the character 
and strength of the inter-particle interaction. In what 
follows we analyze how periodic in r variation of the po- 
tential superimposed on the 1/r 2 decay influences diffu- 
sion. We show that for systems of particles interacting 
via oscillating electron-gas-mediated forces or via forces 
corresponding to the 12 — 6 Lennard- Jones potentials the 
diffusion process depends sensitively on the ratio between 
the distance r at which the potential has a minimum and 
the substrate lattice constant. This ratio allows to dis- 
tinguish between the commensurate and incommensurate 
type of diffusion kinetics. We show that the diffusion co- 
efficient depends periodically on the characteristic length 
of the potentials under investigation. 

The paper is organized as follows. In Section II the 
approach to the diffusion coefficient calculation is shortly 
described. Section III contains description of results for 
electron-gas mediated, oscillation potential , and then for 
12 — 6 Lennard- Jones potential. Section IV summarizes 
main results of the work. 



II. MODEL 

System of N particles interacting via long range forces 
is distributed homogeneously over a one-dimensional sub- 
strate of length L with a lattice constant a. The inter- 
action of two particles at the lattice positions U and lj 
contributes the potential energy s(ali — aL) to the total 
energy of the system. Following Ref. |28||29| we consider 
systems with the pair potential energy e(r) decreasing 
rapidly with r. This justifies to neglecting the next- 
nearest-particle interactions and accounting only for pair 
interactions between neighbors no matter how large the 
intra-pair separation al is. The total interaction energy 
of the system is nie(al) where ni is the number of 
nearest neighbors pairs of length I (in units of a) and 
only Vs satisfying the condition Ini = L are admitted 
in the sum. In a grand canonical ensemble approach we 
let I to vary from to oo and keep, instead, the system 
under fixed external pressure P (in ID it is just an ex- 
ternal force) which is determined by the condition that 
the mean nearest neighbor pair length (I) is equal to the 
inverse of the actual coverage 9 = N/L. In such case a 
probability of a pair of a length I i a 28 i 29 

p,(P,T) = Z^(P,T)e-^ p \ (1) 

where 



and 

oo 

Zx{P,T) =Y,e~ m ' P) (3) 
i=i 

is a single nearest neighbor-pair isothermal-isobaric par- 
tition function. 

Eqs. ([I} through allow to determine the thermody- 
namic properties of the system. In particular, the equa- 
tion of state, relation between the coverage, pressure and 
temperature is obtained by evaluating the mean nearest 
neighbor pair length 

and identifying it with 1/8. In the low temperature limit 
the main contributions to this sum come from one or at 
most two terms only^ 8 -. 

Collective diffusion of the system is modeled by a ki- 
netic lattice gas with the particle hopping rates depend- 
ing on the actual potential energy of the particle. The 
potential energy landscape is build by the static potential 
due to the substrate, as experienced by a single particle, 
and by interactions of the particle with its neighbors. 
Time evolution of the system is controlled by a set of 
master rate equations for the probabilities V({c},t) that 
a microstate {c} of a lattice gas occurs at time t 

= E mm, t) - w{W}, {c})v{{c}, t )\ . 

The microstate {c} is understood as a set of variables 
specifying which particular sites in the lattice are oc- 
cupied and which are not. W({c'},{c}) is a transition 
probability per unit time (transition rate) that the mi- 
crostate {c} changes into {c'} due to a jump of a particle 
from an occupied site to an unoccupied neighboring site. 
Microstates {c} and {c'} differ here only by the position 
of a single particle, the one which jumped. For thermally 
activated jumps the hopping rate depends on the differ- 
ence between the energy of the system when the hopping 
particle is at an intermediate position between the sites 
engaged in the jump and the energy of the system when 
the particle is in its initial position. The only contri- 
butions that do not cancel out in the difference is the 
energy of the hopping particle in its initial position and 
its energy in the activated state at the top of the poten- 
tial barrier which it jumps over. For the particle hopping 
from the adsorption site specified by a pair of integers 
(I, s) (i.e. with the nearest neighbors of adsorbed parti- 
cles being at a distance al and as, respectively, to its left 
and right) to a neighboring site (/', s') — (I ± 1, s =p 1) the 
potential energy at the initial adsorption site is 



e(l,P) = e{al) + aPl 



(2) 



E A = E A + e(al) +e(as), 



(6) 
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where E\ is static potential energy at given site due to 
the interactions with the substrate. The hopping rate 
can be written as 

W({c}, {c'» = Wlf s , = W a e- i A ''°'- e{al) - e{as) ] , (7) 

where W° — v exp[—f3(Eg — E^)] is a hopping rate for 
an isolated (i.e. non-interacting) particle and v is an in- 
trinsic attempt frequency. A;; s s , is the amount by which 
the potential energy E® B of the hopping particle at a 
bridge site between its initial and the final position is 
modified by interactions with the neighbors at each its 
side. We parametrize a microstate {c} as {c} = [X 7 {to}] 
by selecting one particle as a reference particle, denot- 
ing its lattice position as X and specifying positions 
{to} — {mi, rri2, . . . , mjv-i} of all remaining N — 1 parti- 
cles with respect to i1j2&. {to} is referred to as a configu- 
ration. Master equations © are linear in set of probabil- 
ities V(X,{m},t) and so their lattice Fourier transform 
with respect to X can be easily done. The result is that 
^-components of the probabilities, P{ m }{k,t), evolve in 
time independently of each other. The rate equations 
for V{ m }{k 1 t) can be expressed in terms of a fc-space 
microscopic rate matrix M(fc) (with rows and columns 
labelled by {to}) containing the individual hopping rates 
and pha se factors like e ±zka . Details can be found in 
Refs. I25ll28l The collective diffusion coefficient is related 
to that eigenvalue —\n(k) > (termed the diffusive 
eigenvalue) of M(k) which vanishes like k 2 in the limit 
k — > 0. This eigenvalue is then estimated from above in 
a spirit of a variational principle^ 5 , as 



to the "normalization" denominator 



X D (k) 



4>- [-M(fc)] ■ 4> 

J 7 ^ 



-Dk 



(8) 



where the — ► stands for the k — > limit and <fi and <fi are 
variational trial left and right, respectively, eigenvectors 
of M(fc) corresponding to the diffusive eigenvalue. It has 
been show n 25 ' 28 that for a homogeneous substrate the 
{m}-th component of the trial left eigenvector has the 
form 



4>{rn}(k) 



3 = 1 



(9) 



and that c/) {m} (k) = P^ } 4> {m} (k), where P^ n} is the 
probability of the configuration {to} in equilibrium. 

We calculate the diffusion coefficient D as a ratio (the 
k — > limit is implied) 



D 



Xd 

k 2 



M(k) 



N{k)k 2 ' 

of the "expectation value" numerator 

no rep 

{m} ,{m'} 
X W { m'}(k)~4>\ rn} (k) 2 , 



(10) 



{rn} 



(12) 



Detailed balance condition was used to derive Eq. (fTT|) 
so each ({to}, {to'}) term in it accounts for transitions 
between {to} and {to'} in cither direction. Therefore, 
each configuration pair ({to}, {to'}) appears in the sum 
in Eq. (jlip only once [as indicated by the comment "no 
rep" above the sum] in order to avoid double counting. In 
the grand canonical ensemble approach, mentioned ear- 
lier, both and M are functions of P, T, and N. We 
note in passing that TV and M are directly related to the 
diffusion coefficient static (or thermodynamic) and ki- 
netic factor, respectively, which the diffusion coefficient 
is customarily factorized int o 17 ' 18 . The former is con- 
trolled only by the static interactions, determining the 
equilibrium properties of the system, while the latter is 
also sensitive to the dynamic interactions within the ad- 
sorbate and the dynamic interactions with the substrate, 
both controlling the rate of an approach to the thermody- 
namic equilibrium. Certain characteristic features of the 
density dependence of the static factor, often being sig- 
natures of an onset of an organization within the system, 
may or may not be compensated by the features present 
in the kinetic factor, resulting in the density dependent 
diffusion coefficient from which such features may be ab- 
sent. This issue for long range repulsive inter-particle 
interaction was examined in detail in Ref. [28l 

It was shown^ 8 that, for the one dimensional system 
with long range interactions, the denominator (|12[) can 
be expressed as^ 8 - 

(I 2 ) - (I) 2 

J\f(k = 0; P, T, N) = N x '- ^ 1 , (13) 



while the numerator can be written as 



(11) 



M(k;P,T, N) 

OO 

= (kafNj2w^ ltS _ lPl (P^)p s (P,n (14) 



The "no rep" restriction in Eq. pip results in only the 
rates of jumps from the left to right to be explicitly 
present in Eq. (|14p [alternatively, expression mathemati- 
cally equivalent to (I14|) with only the right-to-left jump 
rates explicitly appearing in it can be used]. 

To evaluate Ai(k) using Eq. (TTJJ the potential energy 
correction due to interactions of the activated particle 
is needed. One of the simplest models accounting for 
the activated particle interactions is obtained by realizing 
that the particle hopping from the adsorption site (Z, s) 
to (I + 1, s — 1) surmounts a potential energy barrier at 
a bridge site situated, approximately, at a distance I + h 
and s — i from its nearest left and right adsorbed particle 
neighbor, respectively, and by evaluating the interaction 
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potential energy at the bridge by using e(al) generalized 
to half-integer arguments. Consequently, the potential 
energy correction due to interactions of a particle at the 
bridge site between (I, s) and (Z + 1, s — 1) site is 



A \+l,s-l = £ ( al + fl §) + £ ( aS - a ^)> 



1, 



1, 



(15) 



which, used in Eq. ([7]), leads to the following hopping 
rate for the left-to-right jumps 



l,s TT/0 — /3[e(oi+o^)+e(as— a^)— e(aZ)— e(a5)l 



2+l,s— 1 



Using (fT6|) in IT4|) yields the kinetic factor 



M(k;P,T,N) = {kaf 



[Zi{P,T)Y 



.1=1 



• (16) 



(17) 



in which the definition of e in Eq. ([2]) is used. Similarly 
like in Zi(P,T) in Eq. ([3]) the main contribution to the 
sum over I in Eq. (|17|) comes from one or at most two 
terms for low enough temperatures but, for a given value 
of P, these terms may correspond to different Vs than 
those most significant in Z\. In general, the sum has to 
be evaluated numerically. 



III. RESULTS 

It has been shown that the potential energy of inter- 
actions of adatoms adsorbed on metallic surfaces like 
Cu/Cu(lll) 4 , Co or Co on Cu(lll) 5 , Fe or Co on 
Ag(lll)£, Fe/Cu(lll)2 and Ce/Ag(lll)k vary with the 
inter-atomic distance r as 

f 2sm(S F ) \ 2 sm(2g F r + 2S F ) 

£{r) = ~ tF {-^r- ) — jj^f — (18) 

where e F and q F are, respectively, the Fermi energy and 
the Fermi wave vector of the surface electrons. Shape 
of the potential energy depends also on the Fermi-level 
phase shift 5 F which for many systems is equal to — tt/2. 
This value will be used further in our calculations. The 
inter-particle distance dependence of the interaction en- 
ergy (JTHJ) is plotted in Fig. Q] using a dashed line and com- 
pared with a long range purely repulsive potential energy 
e(r) = a/r 2 used in Refs. l28ll29l . The value of the parame- 
ter a used in Fig.[T]has been chosen in such a way that the 
repulsive potential energy curve forms an upper envelope 
of electron-gas-mediated potential energy (fT5| . Strength 
of both types of interactions decays with the distance like 
1/r 2 , however, whereas the forces corresponding to the 
a/r 2 potential energy are repulsive at any inter atomic 
separation, the potential energy (fig]) oscillates, generat- 
ing attractive forces at some inter-particle distances and 
the repulsive ones at others. Attractive forces in the sys- 
tem lead to the creation of stable structures at the surface 
and should affect the dynamic properties of thesystem. 
We compare here the diffusion kinetics in adsorbates with 



both types of inter-particle interactions and, in addition, 
in systems in which the interactions correspond to the 
12 — 6 Lennard- Jones potential energy 



.^,-^[(2) B -(2)» 



(19) 



typical for interactions between neutral atoms. The pa- 
rameters clj and a decide about depth and position of 
the single potential minimum. The Lennard- Jones po- 
tential has a repulsive oc 1/r 12 wall at short distances, 
much steeper than 1 jr 2 , as it can be seen in the Fig. Q] 
At larger distances this potential is attractive and, de- 
caying like 1/r 6 , it is weak in comparison with any of 
the other two. It will be shown, however, that the exis- 
tence of this attractive part is sufficient for the diffusion 
kinetics in the system with the Lennard- Jones interac- 
tions to be qualitatively similar at some densities to that 
in systems with oscillating interactions. 



A. Diffusion in systems with electron-gas-mediated 
interactions 

It is knowr*2°. that a system of particles with long-range 
unscreened repulsive interactions orders at T — at den- 
sities (coverages) given by a rational fraction (smaller 
than 1) when the interaction potential decays faster than 
1/r. The coverage plotted against external potential has 
a fractal form called a devil staircase^. It is also known^E 
that systems with interactions decaying faster than 1/r 
but additionally screened to nearest neighbors no matter 
how far they are (as described in Sec. [TTJ) also form or- 
dered T — phases at coverages equal to 1 /n where n is 
any natural number. The "phase diagram" (density vs. 
pressure or the chemical potential) in this case is not as 
complicated as that for the unscreened interactions — it 
has no fractal structure. For such a system it is possible 
to calculate all static properties 2 ^ as well as to investigate 
fully the collective diffusion kinetics^ for many models 
of microscopic kinetics. In general, the collective diffu- 
sion coefficient peaks at sufficiently low temperatures at 
densities at which the system orders (except when sharp 
drops in compressibility are compensated fully by kinet- 
ics for very special models of microscopic kinetics 28 ). For 
these densities and for densities around them (at finite 
temperatures) the diffusion is very fast, meaning that 
particles have an ability to rearrange quickly when the 
density changes even by a small amount. 

When the interactions are purely repulsive, like for 
e(r) — a/r 2 , the rearrangement into an ordered phase 
occurs for the entire lattice gas at once, the phase tran- 
sition is of the second order, and formation of ordered 
domains surrounded by disordered regions with lower or 
higher densities is not possible. Such possibility opens for 
systems with interactions like in Eq. (|18p which, depend- 
ing on inter-particle distances alternate between repul- 
sion and attraction. When particles attract each other at 
some distances it is possible to order the system locally 
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FIG. 1: (color online) Inter-particle interaction potential en- 
ergies as a function of the inter-particle distance r. Thick 
solid line (red online): e(r) = a/r 2 with a/a 2 = 0.33; dashed 
line (green online) : the oscillating potential energy in Eq. (JTHJ) 
with qFd — 0.7, = 0.36 and 8f = — tt/2; thin solid line 
(blue online): the 12 — 6 Lennard-Jones potential energy in 
Eq. (TjIJJ with a/a = 2.455 and e L j = 0.031. Energies are 
measured in arbitrary energy units. 



even if globally the density is too low for tha t 11 ! 15 . In 
Fig. [J] we compare the coverage dependence of the dif- 
fusion coefficient plotted at two different temperatures 
for the interaction potential energies plotted in Fig. [1] 
(with the same parameters as there) . It can be seen that 
the character of the curves corresponding to the interac- 
tion potential energies e — a/r 2 and the oscillating one 
[given in Eq. (|18p ] is similar at the same temperatures. 
The only significant qualitative difference between these 
two cases is seen around 9 = 1/3 at the lower of the 
two temperatures considered: the diffusion coefficient for 
the oscillating interaction exhibits a peak almost invisi- 
ble for the system with a purely repulsive interaction at 
the same temperature despite the fact that for purely re- 
pulsive interaction the system does order at 9 = 1/3 at 
T — 0. It indicates that the presence of the inter-particle 
attraction for the system with oscillating interaction al- 
lows for an ordering corresponding to 9 — 1/3 already 
at temperatures much higher than those are needed for 
such ordering with the repulsive long range interactions 
only. 

In fact, r = 3a is the preferred inter-particle distance 
(based solely on on the total interaction energy consid- 
erations) for the potential energy parameters selected in 
this example, so local ordered domains are preferentially 
formed which correspond to a local coverage 9 = 1/3. 
Attraction, creating the potential energy minimum, is a 
necessary condition for the formation of the local ordered 
domains. High diffusion coefficient within such domains 
aids in their formation because collective diffusion of par- 
ticles effectively controls the system ability to create or- 
dered phase o 4 ' 5 ' 11 ! 15 . With purely repulsive interactions 
the ordering at 9 = 1/3 is possible only globally so the 



10000 



1000 















/ 


******* 


- 















0.2 0.4 0.6 0.8 1 

e 



FIG. 2: (color online) Coverage dependence of the collective 
diffusion coefficient at temperatures f3 = 1/ksT = 60 (higher 
curves) and (5 = 20 (lower curves) for the interaction potential 
energies plotted in Fig.[Tj Thick continuous lines (red online): 
for e(r) = a/r 2 ; dashed line (green online): for the oscillating 
potential energy in Eq. (|18[1 : thin solid line (blue online): for 
the 12 — 6 Lennard-Jones potential energy in Eq. (|19p . The 
interaction potential energy parameters are the same as in 
Fig. Q] and 1//3 is expressed in the same arbitrary energy units 
as the interaction potential energy. 



diffusion coefficient peak, very sharp at T = 0, is easily 
smeared out by thermal fluctuations. We can general- 
ize our conclusions here: structures build by attractive 
inter-particle interactions are capable of creating local 
domains of ordered phases. Value of the ratio of the po- 
tential energy spatial oscillation period to the substrate 
lattice constant determines which of the diffusion coef- 
ficient peaks (present for T = at 9 = 1/n for purely 
repulsive long range interactions) is amplified or attenu- 
ated. 

Analyzing Fig. [5] a bit more in detail for the case 
of lower temperatures we see that at coverages below 
9 = 0.5 both the oscillating and the strictly repulsive 
interactions result in the diffusion coefficient of roughly 
the same magnitude except, however, around 9=1/3 
where the broad peak for the oscillating interactions re- 
sults in a somewhat faster diffusion for this case. Both 
types of interactions lead to a sharp increase of the diffu- 
sion coefficient at 9 = 1/2 - another coverage of a T = 
ordered phased. For 9 > 0.5 the diffusion is more effi- 
cient for the purely repulsive interactions than it is for 
the oscillating ones. This can be traced back to the fact 
that the purely repulsive interaction has much stronger 
repulsive core at r — a than the oscillating interaction 
has (c.f. Fig. [I]). All these features are largely washed 
out by thermal fluctuations at the higher temperature of 
the two considered here. 

We show also in Fig. [2] the diffusion coefficient for the 
Lennard-Jones interaction. Here, the interactions favor 
the inter-particle distance of about 3a, almost the same 
as the oscillating interaction does so it is not surpris- 
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FIG. 3: (color online) Inter-particle potential energies (I18[) 
as a function of the inter-particle interaction distance r for 
several values of of (from 0.3/a up to 2.5/a) and e F = 0.36 (in 
the same arbitrary energy units as in the preceding figures.) 

ing that at low concentrations both interactions result 
in almost the same value of the diffusion coefficient. At 
higher concentrations, however, the diffusion kinetics is 
controlled by a very steep repulsive core at short dis- 
tances, much steeper for the Lennard- Jones than for any 
of the remaining two interactions. Consequently, the dif- 
fusion coefficient increases to very high values already at 
9 = 1/3. More features unique to the Lennard-Jones 
interaction will be discussed in Sec. IIII Bl 

In general, the character of the density dependence 
of the collective diffusion coefficient strongly depends on 
how closely the minima and the maxima of the oscillat- 
ing interaction potential energy match the distances be- 
tween particles which occupy the substrate lattice sites. 
In Fig.[3]the interaction potential energies for several val- 
ues of qp are shown. It can be seen that the position of 
the first minimum moves towards higher inter-particlc 
distances with decreasing qp. Geometrically, the mini- 
mum inter-particle distance possible is a, the substrate 
lattice constant, and the overall character of the cover- 
age dependence of the diffusion coefficient is to a major 
extent determined by the character of the interaction at 
this particular inter-particle distance. 

If at the separation a the inter-particle interaction is 
strongly repulsive, as it is for qp — 0.3/a, 0.35/a, 0.7/aor 
1/a, then the diffusion coefficient behaves as a function 
of coverage (for coverages between 9 ~ 1/2 and 1) simi- 
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FIG. 4: (color online) Density dependence of the diffusion 
coefficient in the system with oscillating inter-particle inter- 
action (|18[) for several values of qF (from 0.3/a up to 2.5/a - 
the same as in Fig. |3J), £f = 0.36 and /3 = 60 (the latter two 
parameters are expressed in the same arbitrary energy units 
as in the preceding figures.) 



larly to that corresponding to a purely repulsive interac- 
tion (as seen already in Fig. [5] for qp = 0.7 /a): it raises 
rapidly when the coverage approaches the value 9 = 1/2 
from below, suffers a kink and then it decreases slowly 
with further increase of 9. This behavior was already 
analyzed for a purely repulsive long-range interaction in 
Fig. 9 of Ref. [28|. It can be traced back to the behav- 
ior of the diffusion coefficient static factor, proportional 
to l/fif (i.e. proportional to an inverse of an isothermal 
compressibility), modified by the behavior of the kinetic 
factor, proportional to M.. For these values of qp the 
particles avoid occupying the nearest neighbor sites so 
for S« 1/2 they preferentially occupy every second lat- 
tice site. A substantial additional pressure is needed to 
compress the system above half occupation of the lat- 
tice, the isothermal compressibility is very low, i.e. the 
static factor goes through a sharp maximum at 9 = 1/2. 
The overall shape of the density dependence of the diffu- 
sion coefficient is determined, however, by both the static 
and the kinetic factor. With the hopping rates given in 
Eq. (|16p the high value for interaction potential energy 
at the inter-particle distance a results in a kinetic fac- 
tor which also increases sharply at 9 = 1/2 but does 
not reach a maximum there. Instead, for 9 > 1/2 vari- 
ations of the kinetic factor almost perfectly compensate 
for the variations of the static factor for these densities: 
the kinetic factor continues to increase (albeit, less dra- 
matically than at 9 = 1/2), goes through a broad plateau 
and eventually decreases sharply as 9 approaches the full 
coverage value of 1. The compensation between both fac- 
tors for 1/2 < 9 < 1 results in the diffusion coefficient 
which only slowly decreases as a function of 9 over this 
coverage interval. 

The structure observed in Fig. 0] at coverages below 
9=1/2 (i.e. at 6 « 1/3 and smaller) for q F a < 0.7 
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may be understood in a similar way. In these cases the 
inter-particle interaction is strongly repulsive not only at 
a separation a but also at the next possible one, 2a, so 
at coverages close to 6 « 1/3 the particles preferentially 
occupy every third site with the resulting compressibility 
minimum at that coverage. We can see in Fig. [4] that 
locally, within the coverage interval 0.3 < 9 < 0.47, the 
behavior of the diffusion coefficient for the interactions 
corresponding to qp a < 0.7 is quite reminiscent of that 
within the interval around 9 = 1/2 and above it (for the 
same values of qp and also qp = 1/a): a sharp increase at 
9 = 1/3 is followed by a slow diffusion coefficient decrease 
with increasing 9 (until the singularity at 9 = 1/2 takes 
over). For still lower gp's similar structures are observed 
in Fig. 2] around 9 = 1/6 and even 1/7 but due to their 
overlap they are not very well resolved at the temperature 
selected for the plots. 

Comparing in Fig. [2] the diffusion vs. coverage curves 
evaluated at the same temperatures for the e(r) = a/r 2 
interaction with the ones for the corresponding to it (as 
defined in Fig. |T|) qp = 0.7 /a oscillating interaction (TT5[) 
we see that the presence of inter-particle attraction for 
the latter results in more pronounced structures in D(9). 
The diffusion maximum around 9 = 1/3, very broad and 
almost unnoticeable for the purely repulsive interaction, 
becomes quite obvious and sharp for the corresponding 
to it oscillating interaction. The 6> = 1/3 diffusion coeffi- 
cient peak is clearly visible in Fig. 2] also for qp = 0.35/a 
and 0.3/ a. The diffusion coefficient increase at some ad- 
layer densities (here 9 = 1/3) is due to the attraction 
felt by particles at specific distances. The overall shape, 
however, of the coverage dependence of diffusion curves 
is always a result of several different, often compensat- 
ing each other factors, like a character of interaction at 
different distances (which affect both the static and the 
kinetic factor) , a relative height of the interaction poten- 
tial energy barrier which a hopping particle must over- 
come between adsorption sites (which enters the kinetic 
factor). It is not possible, therefore, to guess from the in- 
teraction potential energy curves alone whether maxima 
of the diffusion coefficient at certain coverages do exist or 
do not. Only analyzing the diffusion coefficient vs. den- 
sity curves we may rationalize the existence of certain 
characteristic features. 

Note that the structure observed around 9 = 1/3 dis- 
appears when qp increases from 0.7 /a to 1/a. Upon fur- 
ther increase of qp also the structure around 9 = 1/2 dis- 
appears. Apparently, the repulsion at the inter-particle 
distance of 2a and a, respectively, ceases to be a factor 
controlling behavior of diffusion. Indeed, we see in Fig. [3] 
that upon increasing qp above qp = 0.7/a the repulsion 
at an inter-particle distance of 2a weakens and is gone for 
qp = 1/a - consequently any structure around 9=1/3 
is no longer present. Then, with further increase of qp 
above qF = 1/a the repulsion at a separation of a also 
weakens and, consequently, no structure around 9=1/2 
is observed for qp = 1.5/ a. It was checked that upon in- 
creasing qp gradually up from 0.7/a the structure around 



9 = 1/3 in Fig. [Devolves initially into the locally concave 
dependence around 9=1/3 and then into the locally con- 
vex one, as seen in Fig. [4] for qF = 1/a. Further increase 
of qp washes out the only remaining structure around 
9 = 1/2, D(9) becomes a structureless concave function 
(as seen for qp = 1.5/ a) and then it evolves into a convex 
one (as seen for qp = 2.5 /a). It is worth noting that the 
evolution of the D(9) from a function with the structure 
at 6 = 1/2 for qp = 1/a through the concave structure- 
less function for qp = 1.5/a and then to the convex one 
for qp = 2.5 /a parallels the evolution of D(9) observed 
in Figs. 3 and 4 of Ref. [23j for the short range interaction 
changing from the strongly repulsive, through the weakly 
repulsive, to the attractive one. Indeed, as seen in Fig. [31 
the inter-particle interaction at a distance a is strongly 
repulsive for qp = 1/a, weakly so for qp = 1.5/a and 
somewhat attractive for qp = 2.5/ a. 

Convex D{9) is characteristic for systems with dom- 
inant attractive inter-particle interactions which cause 
them to bond and to form clusters. When particles stay 
on average at distances that minimize the total interac- 
tion energy, the jump rates are lowered according to ([7]) 
resulting in a decrease of the diffusion coefficient kinetic 
factor. Diffusion for qp = 2.5/ a oscillating interaction 
decreases generally as a function of density. The de- 
crease, quite fast at low coverages, slows down at higher 
ones. 

We have argued above that a sudden change of slope 
of D(9) observed in Fig. [5] at certain characteristic cov- 
erages can be understood as due to strong inter-particle 
repulsion at very short distances. This, however, cannot 
be the only reason for such a behavior because, except 
at 9 = 1/2, such sudden slope changes are much less 
pronounced for purely repulsive interactions, certainly 
the attraction must also play a role. We note also that 
all D(9) curves in Fig. [4] for qp < 1/a start as convex 
functions of coverage already at 9 = and the charac- 
teristic smallest coverage at which they suffer a sudden 
decrease of its slope is 9i = 1/7,1/6,1/3, and 1/2 for 
qFa = 0.3, 0.35, 0.7, and 1, respectively. In each case, the 
corresponding value of a/9i is equal to the inter-particle 
distance at which the potential energy (jTSJ) with the ap- 
propriate value of qp has the first deepest minimum in 
Fig. O The convex character of D{9) at coverages 9 < 9\ 
is a signature of the attractive character of the inter- 
particle interactions at distances beyond the minimum, 
i.e. the distances larger than a/9i. At 6 w 9i the sys- 
tem tends to order with inter-particle distances being, 
on average, equal to a/9i. For 9 > 9i the average inter- 
particle distances are small enough for the repulsion be- 
ing the dominant interaction and ordering is still possi- 
ble at some coverages resulting in D(9) going through 
the second maximum at 9 = 1/3 (for qp > 0.35/a) and 
9=1/2 (for qp > 0.7/a). In short, for a particular oscil- 
lating interaction (i.e. particular qp), the smallest cover- 
age 9i at which D{9) ceases to be a convex function is a 
boundary between the coverages (9 < 9i) for which the 
inter-particle attraction dominates and the ones (6 > 6{) 
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for which the repulsion does dominate. For qp = 2.5/a 
the interaction is effectively equivalent to a short range 
attraction (the first minimum of the interaction potential 
energy is at r < a) while for qp — 1.5 /a the repulsion 
at r = a competes with attraction at r = 2a. For this 
intermediate case the character of D{9) resembles, as ob- 
served earlier, the one appropriate for the weak repulsive 
interactions. 

Character of the inter-particle interaction at closest 
possible separations, i.e. repulsion, attraction, or the po- 
tential energy minimum, determines the overall shape of 
the D{9) dependence. For small ofo's, the particles that 
reside at closest possible separations from their neigh- 
bors experience the interaction induced repulsion. With 
increasing qpa they find themselves first at the interac- 
tion potential energy minimum, then experience attrac- 
tion and then the interaction energy maximum. The cy- 
cle repeats with further increase of qpa. Such a cycle 
of successive repulsions and attractions at a given sep- 
aration should lead to a non-monotonic dependence of 
the diffusion coefficient on the parameter qp at a par- 
ticular fixed coverage: we expect that after a monotonic 
qpa dependence of D for qpa up to about 1 the diffu- 
sion coefficient should pass through a series of minima 
and maxima as qpa further increases. Indeed, this is ob- 
served in Fig.[5]in which D is plotted against qpa for fixed 
coverage. Following the curve for 8=1/3 we observe a 
monotonic decrease of D until a minimum is reached for 
qpa w 1. With further increase of qpa the maximum 
is reached for qpa ~ 1.4, followed by a broad minimum 
around qpa w 2.13, a maximum around qpa 3.55, and 
a minimum again, at qpa w 5. The oscillations of D con- 
tinue with increasing qpa but their amplitude decreases 
mirroring the decreasing amplitude of oscillations of the 
interaction potential energy s(r = const) with increasing 
q F a [c.f. Eq. [fTB"])]. For 9 = 1/2 the qp dependence of 
D is very similar to that for 9 = 1/3 except that the 
former starts with much higher value for small qpa, has 
an inflection point around qpa « 1.4 rather than a nar- 
row minimum followed by a narrow maximum. Beyond 
qpa ks 2 both curves follow closely each other. 
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FIG. 5: (color online) Diffusion coefficient D for — 60 as a 
function of a parameter qF of the oscillating interaction poten- 
tial energy (I18|l with = 0.36 for coverages 6 = 1/?, (dashed 
line) and 8 = 1/2 (solid line). €f and are expressed in 
the same arbitrary energy units as in preceding figures. 







c=2.3a 




a 


=3.4a 










o=2.2a / 

co / 
pi / 


s<* — _____o=1Ja 
















a=1 .0a 
a=0.0 










c 


=0.7a 




a=0.89a 



0.2 0.4 0.6 0.8 1 

e 



FIG. 6: (color online) Coverage dependence of diffusion co- 
efficient for the Lennard-Jones inter-particle interaction (|19|) 
for j3 = 20, €lj = 0.031 and several values of a. €lj and f}~ 
are expressed the same but arbitrary energy units. 



B. Diffusion of particles interacting via 
Lennard— Jones potential 

Neutral particles are known to interact via long-range 
Lennard-Jones potential energy (flT?]) . Its repulsive core 
falls off more rapidly than 1/r 2 (in fact, it falls as 1/r 12 , 
reaches an equlibrium position minimum at r m ; n = 2 1 / 6 cr, 
and is attractive at larger distances, falling off like 1/r 6 ). 
A staircase of ordered phases at coverages 6 = 1/n is 
supported by such an interaction (similarly like for a 
purely repulsive interaction), the most prominent ones 
being 9 = 1/2 and 1/3, and due to huge repulsion at 
short distances the diffusion coefficient of particles or- 
dered in such phases is abnormally large in comparison 
to that when the system is not ordered. We plot in Fig.[S] 



the coverage dependence of the diffusion coefficient at 
fixed temperature in a system with Lennard-Jones inter- 
actions corresponding to several values of a (the other 
parameter in Eq. p^|) . e^j is not varied). The overall 
shape of the presented curves is determined by an inter- 
play between two length parameters in the system: the 
distance r m i n = 2 x ^a below which particles repel each 
other, and the substrate lattice permitted minimum sep- 
aration a between interacting particles. 

We can systematically analyze shapes of the D(9) 
curves in Fig. [5J For a = there are no inter-particle in- 
teractions (except for site blocking) and D(9)/W°a 2 = 1 
as seen in Fig. [5J For a = 0.89a we have r m ; n = a so 
when < a < 0.89a the interaction between particles 
separated by a is attractive. Consequently, it is energeti- 
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cally preferable for particles to be as close to each other as 
possible and for such values of a the D{9) dependence is 
convex (as seen for a = 0.7a and 0.89a in Fig. [6]), typical 
for systems for which the inter-particle attraction dom- 
inates (c.f. Fig. 4 in Ref. H). For a > 0.89a the inter- 
action at the closest separation r — a becomes repulsive 
and, indeed, for a = a, the D(9) dependence becomes 
concave resembling qualitatively D{9) for weakly repul- 
sive short range interaction in Fig. 3 of Ref. 23j. With 
further increase of a the repulsion at r = a increases: for 
(7 = 1.1a we note a characteristic maximum of D around 
6=1/2 also observed already for stronger short range 
interactions in Fig. 3 of Ref. [23. For a = 1.3a the repul- 
sion at r = a is strong enough to induce a preferential 
occupation of every second site (note that the interaction 
at the r = 2a separation is still weakly attractive in this 
case) causing the diffusion coefficient to suddenly raise 
by many orders of magnitude as 9 approaches 1/2 from 
below. 

With further increase of a we reach, at a = 2a/2 1 / 6 = 
1.78a, the point beyond which the interaction at r = 
2a is no longer attractive. For a somewhat larger than 
1.78a we have a very strong repulsion at r = a and a 
much weaker one at r — 2a. We see in Fig. [H] that for 
a = 2.2a a local maximum of D(9) develops for 9 w 1/3 
hinting at a preferential occupation of every third site 
around this coverage due to repulsion for r = 2a, followed 
by a sharp, almost discontinuous raise of D at 9 = 1/2 
due to an extremely strong repulsion at r = a. With 
further increase of a the repulsion at r = 2a becomes 
stronger so the structure around 9 = 1/3 becomes as 
sharp as that around 9 = 1/2. The examples are curves 
for a = 2.3a, 2.6a and 2.8a in Fig. [Sj For the latter, 
however, the interaction at the separation r = 3a is also 
repulsive. In fact, with further increase of a features 
similar to those around 6 — 1/2 and 1/3 develop also 
around 6 = 1/4, as seen for a = 3.4a. This is because for 
a > 3a/2 1 ^ e = 2.67a the interaction at the separation 
r = 3a becomes repulsive and, when strong enough, it 
leads to a preferential occupation of every fourth site for 
9 1/4. This structure is not present yet for a = 2.8a 
at the temperature selected in Fig. [5] 

The character of the coverage dependence of the diffu- 
sion coefficient changes qualitatively every time when the 
interaction between particles becomes repulsive at any 
of the inter-particle separations permitted by the sub- 
strate lattice: i.e. every time when a increases through 
na/2 1 / 6 . For a such that n/2 1 / 6 < a /a < (n + l)/2^ 6 
the particles separated by distances r — a, 2a, . . . ,na 
repel each other (with the repulsion being stronger for 
shorter separations) while they attract each other for 
separations r = (n + l)a and larger. In such case one 
expects at T « a sharp increase of the diffusion co- 
efficient at coverages 9 = 1/2, 1/3, . . . , l/(n + 1) due to 
the low temperature structural organization of the adsor- 
bate at these coverages to minimize the total interaction 
energy in the system. In our runs done at lowest temper- 
atures for which the calculations are feasible we see such 
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FIG. 7: (color online) Dependence of the diffusion coefficient 
at 6 = 1/5 on the Lennard- Jones interaction parameter a for 
(3 = 20 and e L j = 0.031 (same as in Pig. EJ. 

structures for 9 = 1/2, 1/3 and 1/4. Still, one expects 
that as the parameter a is varied one should observe a 
non monotonic oscillatory changes of the diffusion coef- 
ficient at low enough coverages. Indeed, this is seen in 
Fig. [7] in which D(9 = 1/5) is plotted as a function of 
a /a. The maxima are noted for a /a approximately half 
way between a = na/2 1 ^ 6 for n = 1,2,3,4,..., i.e. for 
cr/a « 1.3,2.2,3.1,..., for which the interaction is al- 
ready strongly repulsive for all inter-particle separation 
up to a, 2a, 3a, . . ., respectively. 

IV. CONCLUSIONS 

Variational approach has been applied to examine col- 
lective diffusion in a one dimensional lattice gas sys- 
tems with two type of long-range inter-particle inter- 
action: the electron-gas-mediated interaction described 
by the oscillating FriedeHike potential energy (|18p. and 
the Lennard- Jones interaction corresponding to the po- 
tential energy (|T9|l . We have discussed the features of 
the coverage (adsorbate density) dependence of the diffu- 
sion coefficient for both these interactions and compared 
them with those investigated in detail for purely repulsive 
long-range interaction^ corresponding to the potential 
energy oc 1/r 2 as well as with the behavior typical for 
short-range repulsive and attractive interactions^. 

In general, at densities above half coverage (9 > 1/2), 
at which the inter-particle repulsion at short distances 
plays a main role, the diffusion coefficient for the repul- 
sive and the oscillating interaction behaves, as a function 
of coverage, qualitatively similarly: the diffusion coef- 
ficient is much higher than without interactions and de- 
pends on coverage relatively weakly. This is true even for 
the Lennard-Jones interaction, except that a very steep 
repulsive core in this case, making creation of a high den- 
sity adsorbate energetically very costly, results in huge 
values of the diffusion coefficient at such densities. 
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At lower densities, the behavior of the coverage depen- 
dence of diffusion coefficient for the Lennard-Jones inter- 
action, repulsive for inter-particle separations r < r m ; n 
and attractive for r > r m - m , is quite easy to under- 
stand. At sufficiently low temperatures, it experiences 
for increasing 9 a finite series of progressively sharper in- 
creases at "critical" coverages d — 1/n (n — n max , n max — 
1, . . . , 2). The first one, for 9 = l/ra max corresponds to 
the largest inter-particle distance r = n max a < r m - m at 
which the interaction between the particles is still re- 
pulsive (i.e. at the separation r = (ri max + l)a > r m ; n 
the interaction is already attractive). As temperature 
increases, the structures at lower densities of this series 
are usually smoothed out. In fact, we observe for the 
Lennard-Jones interaction a delicate interplay between 
two length scales: the lattice constant a which deter- 
mines what actual distances between particles are possi- 
ble, and the characteristic length of the interaction po- 
tential, r m j n , separating the short range repulsion from 
the attraction at larger distances. When r min < a then 
the repulsive core of the interaction is irrelevant and the 
diffusion coefficient exhibits features similar to those ob- 
served for short-range attractive interactions. 

Features of diffusion for the oscillating interaction are 
somewhat more difficult to explain in simple terms due 
to the existence of multiple interaction potential energy 
minima and, consequently, alternating regions of inter- 
particle attraction and repulsion. Still, a sudden raise 
of the diffusion coefficient at 9 = 1/2 has origin similar 
to that for such structure for the purely repulsive and 
Lennard-Jones interaction. The structures for 9 < 1/2 



are, however, masked by the influence which attraction 
alternating with repulsive interactions has on diffusion. 

We must note also that interpreting features of the 
diffusion coefficient in terms of the features of the inter- 
particle interaction is deficient in that respect that it nec- 
essarily is limited to the interpretation of the kinetic phe- 
nomenon (diffusion) in terms of the static properties of 
the system, i.e. in terms of the static (or thermodynamic) 
factor commonly used in theories of diffusion. The ki- 
netic factor is known, however, to be capable of compen- 
sating for very often drastic behavior of the static factor 
as a function of density^. This is also the case here: 
a slow variation of the diffusion coefficient as a func- 
tion of coverage for several curves in Figs. 2] and [5] for 
9 immediately larger than 1/2, 1/3, and 1/4 is the result 
of a delicate compensation between a sharp drop of the 
static factor (which at these particular coverages suffers a 
sharp cusp-like maximum corresponding to a sharp drop 
in an isothermal compressibility) and a strong increase of 
the kinetic factor when the coverage is increased through 
these values. 
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